function  [f, df] = dace_regpoly2(S)
%Second order polynomial regression function
% Call:    f = dace_regpoly2(S)
%          [f, df] = dace_regpoly2(S)
%
% S : m*n matrix with design sites
% f =  [1 S S(:,1)*S S(:,2)S(:,2:n) ... S(:,n)^2]
% df : Jacobian at the first point (first row in S)

% hbn@imm.dtu.dk
% Last update September 4, 2002

[m n] = size(S);
nn = (n+1)*(n+2)/2;  % Number of columns in f
% Compute  f
f = [ones(m,1) S zeros(m,nn-n-1)];
j = n+1;   q = n;
for  k = 1 : n
    f(:,j+(1:q)) = repmat(S(:,k),1,q) .* S(:,k:n);
    j = j+q;   q = q-1;
end

if  nargout > 1
    df = [zeros(n,1)  eye(n)  zeros(n,nn-n-1)];
    j = n+1;   q = n;
    for  k = 1 : n
        df(k,j+(1:q)) = [2*S(1,k) S(1,k+1:n)];
        for i = 1 : n-k
            df(k+i,j+1+i) = S(1,k);
        end
        j = j+q;   q = q-1;
    end

end

return
